We estimate infection fatality rates and hospitalization rates from [Verity et al. 2020] from data from the Wuhan outbreak. These are corrected for both asymptomatic infection and underreporting of symptomatic cases. We assume that between 20 - 40 % of people are cummulatively infected.
Figure 1. Predicted IFR by age group
Figure 2. Predicted hospitalization rates by age group.
Importantly the IFRs are estimated from the Wuhan context and could be thought of as a lower bound. Thus, the hospitalization rates it would be important to know how many people die in the absence of ventilators, etc.
---
title: "Estimating the burden of COVID-19 in African countries"
output:
flexdashboard::flex_dashboard:
vertical_layout: fill
orientation: row
source_code: embed
theme: simplex
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(message = FALSE, warning = FALSE)
# Load Packages
library(htmltools)
library(DT)
library(dplyr)
library(data.table)
library(leaflet)
library(rgdal)
library(ggplot2)
library(plotly)
# set parameters here
p_infected <- 0.2 # 20% cummulative infections
age_brackets <- c("A0004", "A0509", "A1014", "A1519", "A2024", "A2529", "A3034", "A3539", "A4044",
"A4549", "A5054", "A5559", "A6064", "A65PL") # names should match colnames!
# hospitalization rates from Verity et al. 2020 table 3
# 10 yr age brackets, hence rep for 5 yr age brackets & % hence divide by 100
hosp_mean <- rep(c(0, 0.0408, 1.04, 3.43, 4.25, 8.16, 11.8)/100, each = 2)
hosp_lower <- rep(c(0, 0.0243, 0.622, 2.04, 2.53, 4.86, 7.01)/100, each = 2)
hosp_upper <- rep(c(0, 0.0832, 2.13, 7.00, 8.68, 16.7, 24.0)/100, each = 2)
# ifr for these brackets from Verity et al. 2020 table 1
ifr_mean <- rep(c(0.00161, 0.00695, 0.0309, 0.0844, 0.161, 0.595, 1.93)/100, each = 2)
ifr_lower <- rep(c(0.000185, 0.00149, 0.0138, 0.0408, 0.0764, 0.344, 1.11)/100, each = 2)
ifr_upper <- rep(c(0.0249, 0.0502, 0.0923, 0.185, 0.323, 1.28, 3.89)/100, each = 2)
# set the names
names(hosp_mean) <- names(hosp_lower) <- names(hosp_upper) <- names(ifr_mean) <- names(ifr_lower) <- names(ifr_upper) <- age_brackets
```
```{r ests, include = FALSE}
# read in data
country_dt <- fread("output/country_dt.csv")
admin_dt <- fread("output/admin1_dt.csv")
# apply to data @ country/admin level
admin_dt[, c("deaths_mean", "deaths_upper", "deaths_lower") :=
.(rowSums(admin_dt[, Map("*", .SD, ifr_mean*p_infected), .SDcols = names(ifr_mean)],
na.rm = TRUE),
rowSums(admin_dt[, Map("*", .SD, ifr_upper*p_infected), .SDcols = names(ifr_mean)],
na.rm = TRUE),
rowSums(admin_dt[, Map("*", .SD, ifr_lower*p_infected), .SDcols = names(ifr_mean)],
na.rm = TRUE)), .SDcols = names(ifr_mean)]
admin_dt[, c("hosps_mean", "hosps_upper", "hosps_lower") :=
.(rowSums(admin_dt[, Map("*", .SD, hosp_mean*p_infected), .SDcols = names(hosp_mean)],
na.rm = TRUE),
rowSums(admin_dt[, Map("*", .SD, hosp_upper*p_infected), .SDcols = names(hosp_mean)],
na.rm = TRUE),
rowSums(admin_dt[, Map("*", .SD, hosp_lower*p_infected), .SDcols = names(hosp_mean)],
na.rm = TRUE))]
admin_dt$prop_ov65 <- admin_dt$A65PL/admin_dt$pop
admin_dt$inc_per100k <- admin_dt$deaths_mean/admin_dt$pop*1e5
# country level
# apply to data @ country/admin level
country_dt[, c("deaths_mean", "deaths_upper", "deaths_lower") :=
.(rowSums(country_dt[, Map("*", .SD, ifr_mean*p_infected), .SDcols = names(ifr_mean)],
na.rm = TRUE),
rowSums(country_dt[, Map("*", .SD, ifr_upper*p_infected), .SDcols = names(ifr_mean)],
na.rm = TRUE),
rowSums(country_dt[, Map("*", .SD, ifr_lower*p_infected), .SDcols = names(ifr_mean)],
na.rm = TRUE)), .SDcols = names(ifr_mean)]
country_dt[, c("hosps_mean", "hosps_upper", "hosps_lower") :=
.(rowSums(country_dt[, Map("*", .SD, hosp_mean*p_infected), .SDcols = names(hosp_mean)],
na.rm = TRUE),
rowSums(country_dt[, Map("*", .SD, hosp_upper*p_infected), .SDcols = names(hosp_mean)],
na.rm = TRUE),
rowSums(country_dt[, Map("*", .SD, hosp_lower*p_infected), .SDcols = names(hosp_mean)],
na.rm = TRUE))]
country_dt$prop_ov65 <- country_dt$A65PL/country_dt$pop
country_dt$inc_per100k <- country_dt$deaths_mean/country_dt$pop*1e5
# From http://leafletjs.com/examples/choropleth/us-states.js
admin1 <- readOGR("output/shapefiles/admin1.shp")
admin1$id_match <- 1:nrow(admin1@data)
country <- readOGR("output/shapefiles/country.shp")
# merge
admin1@data <- left_join(admin1@data,
select(admin_dt, starts_with("death"), starts_with("hosp"), starts_with("pseek"),
pop, admin_name, id_match, admin_type, prop_ov65,
inc_per100k))
```
Sidebar {.sidebar}
======================================================================
The goal of this project is to map the potential burden of COVID-19 in African countries. Currently we are focusing on demography, but hope to incorporate other factors such as healthcare capacity.
A summary of our approach is described in the __Methods__ tab. We use demographic data from [WorldPop](https://www.worldpop.org/geodata/summary?id=1276), administrative data from [The Malaria Atlas Project](https://malariaatlas.org) accessed through the R package [`malariaAtlas`](https://cran.r-project.org/web/packages/malariaAtlas/index.html). We also have summarized data from [Alegana et al] on probability of seeking treatment and travel times to health facilities. Download these data and the associated metadata here. We take data from [Verity et al. 2020](https://www.thelancet.com/journals/laninf/article/PIIS1473-3099(20)30243-7/fulltext#seccestitle230) & also assume cummulative infection rate
To explore these assumptions, check out our shiny app. [here](https://marjoleinbruijning.shinyapps.io/covid19-burden-africa/).
This project is an offshoot of this work by [Miller et al](https://github.com/ianfmiller/covid19-burden-mapping) mapping burden and hospital capacity in the US.
All code & data are available [here].
Mapping burden
======================================================================
Row {data-height=1000}
-----------------------------------------------------------------------
### Burden mapped to country & Admin 2 level (where available)
```{r}
bins <- c(0, 20, 40, 60, 80, 100, 150)
pal <- colorBin("Reds", domain = admin1$inc_per100k)
labels <- sprintf(
"Country: %s
Admin name (type): %s (%s)
Pop: %s
% pop over 65: %0.2g
Estimated deaths: %0.2f (%0.2f - %0.2f)
Estimated hospitalizations: %0.2f (%0.2f - %0.2f)
Estimated reporting of fevers to hospitals: %0.3g",
admin1$name_0, admin1$admin_name, admin1$admin_type,
format(admin1$pop, big.mark = ",", scientific = FALSE, digits = 0), admin1$prop_ov65*100,
admin1$deaths_mean, admin1$deaths_lower, admin1$deaths_upper,
admin1$hosps_mean, admin1$hosps_lower, admin1$hosps_upper,
admin1$pseektrthosp10x10
) %>% lapply(htmltools::HTML)
leaflet() %>%
setMaxBounds(-25.35875, -40.37063, 63.5003, 37.54327) %>% # from bbox(admin2)
addProviderTiles('CartoDB.Positron') %>%
addPolygons(data = admin1,
color = "black", weight = 0.001, smooth = 0.3,
fillColor = ~pal(inc_per100k),
fillOpacity = 0.7,
dashArray = NULL,
highlight = highlightOptions(
weight = 3,
color = "red",
dashArray = NULL,
fillOpacity = 0.7,
bringToFront = TRUE),
label = labels,
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto")) %>%
addLegend(pal = pal,
values = admin1$inc_per100k, opacity = 0.7,
position = "bottomright", title = "Deaths per 100,000 persons") %>%
addPolygons(data = country, color = "#444444", weight = 2, fill = FALSE)
```
Row {data-height=600}
-----------------------------------------------------------------------
### Relationship between admin/country level ests of burden
```{r}
country@data <- left_join(country@data,
select(country_dt, iso, starts_with("death"), starts_with("hosp"),
starts_with("pseek"), pop, prop_ov65))
p <- ggplot() +
geom_point(data = filter(country@data, !is.na(hosps_mean)),
aes(x = reorder(name_0, prop_ov65), y = hosps_mean/pop*1e5,
color = prop_ov65), shape = 15) +
geom_point(data = filter(admin1@data, !is.na(hosps_mean)),
aes(x = name_0, y = hosps_mean/pop*1e5, color = prop_ov65),
shape = 22, alpha = 0.5, size = 0.5) +
labs(y = "Hospitalizations per \n 100,000 persons", x = "") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_color_distiller(palette = "PuRd", name = "Proportion of pop \n over 65")
ggplotly(p)
```
Methods
======================================================================
We estimate infection fatality rates and hospitalization rates from [Verity et al. 2020] from data from the Wuhan outbreak. These are corrected for both asymptomatic infection and underreporting of symptomatic cases. We assume that between 20 - 40 % of people are cummulatively infected.
Figure 1. Predicted IFR by age group
Figure 2. Predicted hospitalization rates by age group.
Importantly the IFRs are estimated from the Wuhan context and could be thought of as a lower bound. Thus, the hospitalization rates it would be important to know how many people die in the absence of ventilators, etc.